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We report a theoretical scheme that enables the calculation of maximally localized Wannier func- 
tions in the formalism of projector-augmented-waves (PAW) which also includes the ultrasoft- 
pseudopotential (USPP) approach. We give a description of the basic underlying formalism and 
explicitly write all the required matrix elements from the common ingredients of the PAW/USPP 
theory. We report an implementation of the method in a form suitable to accept the input electronic 
structure from USPP plane-wave DFT simulations. We apply the method to the calculation of Wan- 
nier functions, dipole moments and spontaneous polarizations in a range of test cases. Comparison 
with norm-conserving pseudopotentials is reported as a benchmark. 
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I. INTRODUCTION 

Wannier functions (WFs) can be obtained by a unitary 
transformation of the extended wavefunctions of a peri- 
odic systemic So far, the effective diffusion/apphcation 
of WFs in electronic structure calculations was hindered 
by their intrinsic non-uniquenessi^ In 1997 Marzari and 
Vanderbilt proposed a useful approach^ii to overcome 
this drawback: the proposed methodology allows one to 
extract, from a selected manifold of bands, the set of 
WFs with the maximum spatial localization, i.e. the 
maximally localized Wannier functions (MLWFs). On 
one hand, MLWFs are attractive because they consti- 
tute a complete and orthonormal basis set in the real 
space. On the other hand, with respect to other nu- 
merical real-space basis sets, they carry also the phys- 
ical information of the starting Bloch functions. In- 
deed, MLWFs may yield the chemical view of molec- 
ular bond orbitals, and they can be exploited for the 
computation of the spontaneous polarization in periodic 
systemSji^^^ becoming very popular to tackle these issues 
in advanced materials iSiSiiS In addition, MLWFs have 
most recently been proposed to calculate the transport 

properties of nano-size conductors connected to external 
electrodes,"d^ii3,i4,i5 

The calculation of MLWFs was originally imple- 
mented"^ for a projection from Bloch orbitals expanded 
on a plane-wave basis set in Density Functional Theory 
(DFT) calculations, within the norm-conserving pseu- 
dopotential (NCPP) framework. Norm-conserving pseu- 
dopotentials, "'^^•-'^^ that allow to neglect core electrons in 
the evaluation of physical observables, are usually charac- 
terized by a high transferability of an element to a variety 
of chemical environments. However, the norm conserva- 
tion in the core region is a strong constraint, that affects 
the computational effort of the DFT calculations. As a 
consequence, only a part of the periodic table elements re- 
sults numerically accessible: some chemical species, such 
as first-row elements (e.g. C, N, O, F) and especially 
transition metals (e.g. Mn, Fe, Co, Ni, Cu) and rare 



earths (e.g. La, Gd, Yb), would require an extremely 
high number of basis functions {e.g. plane waves), in 
order to be described with a satisfactory accuracy. Un- 
fortunally, the typical systems of interest in nanoscience 
and specifically in molecular electronics contain atoms of 
those critical species. 

This problem was brilliantly solved within a frozen- 
core approach by introducing the ultrasoft pseudopo- 
tentials (USPPs)fi^ that relax the norm conservation 
constraint, by compensating with a pseudized charge. 
Since the computational requirement in the calculation 
of the DFT Bloch functions affects also the evaluation of 
the MLWFs, the generalization of the original Marzari- 
Vanderbilt's procedure to the case of USPPs becomes 
necessary, in order to make the MLWFs a powerful tool 
to tackle realistic systems of high technological and fun- 
damental relevance. This is the focal aim of the paper. 

In pursuing this extension, we also noted that recent 
studies mapped^^-^° the USPP procedure into the PAW 
theory. The latter^^-^^ approach was developed by Blochl 
to combine an all-electron description of the system with 
the simplicity of the frozen-core pseudopotential meth- 
ods. The reader is referred to the original articles for a 
description of the PAW scheme and for the matching be- 
tween PAW and USPP. The USPP can indeed be recast 
in the PAW formalism as an approximation, as we will 
discuss in the following. Given this equivalence, we de- 
veloped a theoretical framework within the PAW theory 
to compute WFs from USPP Bloch orbitals, as we show 
by expressing the necessary matrix elements. 

The paper is organized as follows: in section ^ we 
first write down the salient quantities of the USPP and 
PAW methods that enter the computation of the WFs, 
and then obtain and discuss the WF computation; in 
section Unl we report the results of the application of our 
method to several test cases, that explore both the chem- 
ical bonding and the electrical polarization in the perti- 
nent cases; finally we draw our conclusions in section llVl 
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II. FORMALISM 



B. PAW and ultra-soft pseudopotentials 



A. Maximally localized Wannier functions 

In this section we give a brief introduction to the theory 
of maximally localized Wannier functions. A more de- 
tailed description can be found in the original papers^*^ 
In the case of an isolated band, Wannier functions can 
be definedii^ as a combination of the Bloch orbitals j-i/'k) 
corresponding to different k— points as follows: 
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where e"^'' is a k-dependent phase factor. This definition 
has been generalized^ to a group of bands leading to the 
expression: 
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Here the extra degrees of freedom related to the phases of 
the Bloch eigenstates are collected in the unitary matrix 
U^. In the one-band case it has been demonstrated that 
a suitable choice of the phases e"^'' leads to WFs which 
are real and exponentially decaying in a real space rep- 
resentation!^ In the many-band caseSi this theorem does 
not hold anymore but the arbitrariness in the unitary 
(gauge) transformation can be exploited. Following 
Marzari and Vanderbilt^ we define a spread functional 
ft [U^] , which gives a measure of the degree of localiza- 
tion of the WF set. It reads: 



(3) 



where (■)„ is the expectation value of a given operator 
on the n-th WF calculated using the gauge transfor- 
mation. It is therefore possible to define the maximally 
localized Wannier functions (MLWFs) as the WFs result- 
ing from Eq. Q by means of the unitary transformation 
that minimizes the spread functional. 
According to the formal analysis of the (f)„ and {f'^)n 
terms given by Blount, ■^'^'^ it is possible to demonstrate 
that the dependence of on the gauge transformation is 
determined only by the so-called overlap integrals M^'^: 
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drul (r) Uk+b,«(r), 



(4) 



Mk,m(r) being the periodic part of the Bloch states 
4'k,m{'r) = e'*' '' Uk,jn(r). The detailed form of the po- 
sition expectation values and of the spread functional 
(including its gradient wrt to U^) in terms of the overlap 
integrals is reported in Appendix^ Since the represen- 
tation of Bloch eigenstates enters only the calculation of 
the M^'^ integrals, these quantities are the main objects 
to deal with when using a USPP or PAW formalism. The 
detailed treatment is reported in section III CI 



The PAW formalism has been introduced by 
BlochlSiiSSiS^ and it has also been demonstrated^ that 
the Vanderbilt ultrasoft pseudopotential (USPP) the- 
Qj,yi8,i9^26,27^28 ^^^^ j-,g obtained within the PAW ap- 
proach. 

Blochl's starting point is to partition the volume of 
the system by setting spherical regions (atomic spheres) 
around each atom. In each sphere two complete sets of 
wavelets^ ({It/)"*^)} and {If/if*)} ) localized in the sphere 
are defined. While the former, when truncated, is in- 
tended to work with all-electron (AE) functions, the lat- 
ter should be smoother and easily representable in plane 
waves. It is therefore possible to introduce a well de- 
fined linear operator mapping one-to-one AE-like func- 
tions into smoother pseudo (PS) functions and viceversa: 



where 



(5) 
(6) 



The index / runs over different atoms {e.g. different 
atomic spheres), while {Pn \ are the projectors-*^ related 
to the pseudo wavelets \4>^i). The T operator acts on the 
pseudized functions to reconstruclii^ the AE ones. 

Matrix elements and expectation values of a generic 
operator A on the physical AE states can be written as: 



One of the main results of Ref. |23| is the explicit expres- 
sion for Ap**, which we report here for the case of local 
and semilocal operators: 



AP' = A + A'''^\ 



(8) 
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The second term (A""^) in the rhs of Eq. © takes into 
account the corrections due to the use of the pseudo func- 
tions instead of the AE ones, and from here on it will be 
called augmentation term. At this point it is useful to de- 
fine the quantities (r) and ql^ (augmentation densities 
and charges respectively) as: 

Ql^ (r) = rit * (r) (r) - C * ■ (r) , (9) 



I dvQl^(v) 



(10) 



Within these definitions, Eq. © for local operators A{y) 
can be recast in a more convenient form: 



drQf,(r)A(r) 
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Setting A to the identity in Eqs. (|8I11() . scalar products 
are given by ('0m | S IV'n*) where the number operator— S 
(that characterizes also the USPP formalismiS) is given 
by: 

5 = ftf = 1 + ^ |/3,,)g(,.(/3,^.|. (12) 

In the same way, by setting A{r') — eS{r' — r) we obtain 
an expression for the density: 

n{r) = nP%r) + n'^'^s {r) , (13) 

mk I ,ij 

where n^''(r) is the density contribution of the pseudo 
wavefunctions. Since we are able to express all the quan- 
tities of interest in terms of the soft pseudo-states, the 
quantum problem can be solved directly in this represen- 
tation. In order to do this, it is necessary to augment the 
Hamiltonian operator: the procedure leads to additional 
terms which have exactly the same role as that of the 
pseudopotentials in standard PW calculations. 

Moving to the USPP framework, the generalization in- 
troduced by Vandcrbilt— is twofold, (i) More than one 
projector per angular momentum channel can be taken 
into account: the inclusion of multiple projectors per 
channel enlarges the energy rangeiSiSi over which loga- 
rithmic derivatives are comparable with the full potential 
case, thus increasing the overall portability of the pseu- 
dopotcntial. (ii) By relaxing the norm-conservation con- 
straint of the pseudo reference-states, the pseudo wave- 
functions are smoothened: thus, the required cutoff en- 
ergy for PW representation can be drastically lowered. 
The fact that properties (i) and (ii) are verified for the 
PAW wavelets (pi establishes the connection between the 
PAW and USPP methods. In fact, (i) is naturally valid 
for the wavelets (/)^^, because these form a basis set, and 
therefore have in principle an infinite number of states 
for each angular momentum channel, (ii) is valid as well, 
and the possibility of non norm-conservation in passing 
from \4>i'^) to l^f*) is accounted by non null qfj terms 
in Eqs. (|9I1U|) . Consequently, the PAW theory for wave- 
function reconstruction can be basically adopted also in 
the case of USPP,^2iS 

While the PAW method is in principle an exact AE 
(frozen-core) approachfS the USPP method adopts a 
further approximation22iS represented by the require- 
ment of pseudizing the augmentation densities Qfj(r) 
[Eq. 0]. Since these terms contain the AE reference 
states, they are not simply writable on a PW basis: the 
USPP pseudization is done to make them suitable for a 
PW representation. This also means that the total den- 
sity from Eq. H13|) would be PW representable within 
such an approach. Even though the augmentation den- 
sities are pseudized, they must capture some features of 
the physical AE density. Consequently, their PW cut- 
off energy may be larger than the one associated to the 
pseudo wavefunction density [^^''(r) in Eq. (|13|l ]. 



C. Maximally localized Wannier functions within 
PAW/USPP 

As mentioned in Sec. Ill Al Bloch wavefunctions enter 
the calculation of MLWF's only through the overlap ma- 
trix elements M^'^, defined in Eq. Q. Therefore, the 
reconstruction of these integrals from the knowledge of 
pseudo (ultrasoft) functions completely solves the prob- 
lem of computing MLWF's within a PAW/USPP formal- 
ism. Once overlap matrices have been calculated, we do 
not longer distinguish whether the parent Bloch wave- 
functions were pseudized or not. 

Being the overlaps [Eq. the matrix elements 

of the local operator e"*''''", we built up the corresponding 
augmentation using Eq. (|11() . Overlaps can be written as: 

= KJ<+^J + (14) 

where we have defined Qjj (b) — J dr Qj^ (r) e~^^'^ and 
/3jj) are k-point symmetrized projectors (Bloch sums). 
Summation over ions / is done in a single unit cell. De- 
tails about the calculation of these quantities are re- 
ported in Appendix IBl We stress that the scalar product 
of lit^"*) functions cannot be simply augmented by the S 
[Eq. H12I) ] operator as those involving |'0^*)'s. In order 
to work with the periodic part of the Bloch func- 

tions, first we have to reconstruct the AE Bloch states by 
means of T and then we can obtain the required |w^',„) 
states by applying the local operator e~^^'^ . Since this 
last operator does not commute with T, we are not al- 
lowed to directly work on [m^"*) with the reconstruction 
operator. In the first scalar product of Eq. (|14|l the num- 
ber operator S — T'^T has been therefore substituted by 
the augmented operator M — e"'*' '' T. We thus need 
to introduce the Fourier transform of the augmentation 
densities Qf^ (b) instead of the augmentation charges ql^ . 
We also note that in the thermodynamic limit the It- 
point grid becomes a continuum, so that b ^ 0. This 
limit leads to identical S and M operators. Therefore, 
within discrete k-meshes the use of S instead of M is an 
approximation expected to give the best performance in 
the limit of a large number of k-points. We will refer to 
this approximation as the thermodynamyc limit approx- 
imation (TLA). We will comment more on its numerical 
aspects in Sec. IIIII 

The problem of calculating Wannier functions within 
USPP has also been faced elsewhere in the litera- 
^yj.g|34j35j36 a first attempt, Vanderbilt and King- 
SmithM extended the calculation of the spontaneous po- 
larization through the Berry phas o^i^i"^ to the USPP pro- 
cedure. Since the Berry phase*" is directly related to over- 
lap integrals, an expression for its calculation is reported 
in Eq. (23) of Ref. 34]. This expression adopts the S 
number operator instead of M. Therefore, the result 
does not completely agree with the one presented here 
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[Eq. (|14(l ]. but it can be read as an approximating for- 
mula having the right thermodynamic Hmit in view of 
the above discussion. Bernasconi and Madde n?^'?^ de- 
rived instead the formahsm for MLWFs with USPP in a 
simphfied approach"^^ vahd only in the case of F-sampled 
supercells. Although the basic ingredients (e. (7. overlaps) 
are the same, our treatment is valid for generic periodic 
systems, and recovers the F-only calculation as a special 
case. Furthermore, the proof given in Ref. is not 
based on a PAW reconstruction as we did. The authors 
generalized the augmentation operator for the density 
[Eq. H13|) ] also to the case of the density matrix and then 
derived the augmentation for the e"**^' '' operator. Here 
instead the are the generators of the reciprocal lattice. 
We note that the result by Bernasconi and Madden is 
coherent with our treatment, while the first one by Van- 
derbilt and King-Smith^^ is not. Following Bernasconi 
and Madden, Thygesen and coworkers'^^ recently derived 
an expression for overlap matrices within USPP for pe- 
riodic systems, considering a F-only supercell containing 
the whole crystal. 

Closing this section we wish to underline that matrix 
elements of the form (V'lcml e~'^'i+*^^ '' |?Ak+q,n) enter also 
other physical problems. A particularly appealing case 
is the calculation of one-particle Green function in the 
GW approximationi^24£ Current implementationgiii^Sii^ 
of the method experience DFT wavefunctions mainly 
through the above defined matrix elements (evaluation of 
the polarizability), while no direct access to the density 
is required. The extension of GW calculations to the case 
of PAW^^^ or USPP is therefore feasible along the same 
lines we presented here for Wannier functions. We note, 
however, that the numerical cost is extremely higher due 
to the larger number of overlaps to be computed. 



\(3f) and Qij{h) are reported in Appendix [51 Since no 
reference to the charge is made, only the very smooth 
wavefunction grid is used throughout the calculation. 
While a linear scaling with the number of plane waves is 
exploited in the first term (pseudo overlaps) of Eq. H14|l . 
scalar products between projectors and pseudo states in 
the second term (augmentation overlaps) are the price 
to pay for introducing USPPs. When we consider that 
the number of /3-projectors Np is of the same order of 
the number of bands Nh (but usually larger by a factor 
between one and two) we see that both pseudo and aug- 
mentation overlap terms have the same scaling, namely 
X A^k X Np\Y- However, the pseudo overlaps turn out 
to have a larger prefactoi4£ and represent the leading 
term. Usually, USPPs allow for a reduction of the PW 
cutoff by a factor of 2 to 3 for first-row elements up to 5 or 
even more for atoms with d- or /-states. This leads to a 
reduction of the PW number Npw by a factor of around 
3 to 10 or more. Even if the scaling wrt A^pw is linear, it 
more than compensates the effort for augmenting over- 
laps and makes the introduction of USPPs numerically 
advantageous. Our experience shows that USPPs avoid 
the creation of bottlenecks in the computation of over- 
laps and make the WF localization the leading part of 
the calculation. 

Finally, we note that in order to give a guess for the 
iterative minimizations involved in the MLWF method, 
it is sometimes required to compute the projections of 
Bloch states onto some starting localized functionsi^ 
The augmentation of scalar products is performed as 
usual accounting for the S number operator [Eq. IT^ . 
Projections on the /3-functions are required as well, but 
they have already been computed and it turned out that 
scaling is linear wrt the PW number as before. Therefore, 
no implications on the above discussion arise. 



D. Numerical details 

In this section we would like to describe some issues 
related to the numerical performance of the method. We 
implemented this formalism in the freely-available WanT 
code,*^ for the calculation of electronic and transport 
properties with WFs. We also took advantage of the 
complete integration of the WanT code with the PWSCF 
packager^ which explicitly treats the DFT problem using 
USPPs. From here on we focus on the USPP formalism, 
even though a large part of the discussion is still valid 
also in the PAW case. The most important advantage 
of using the USPP construction for Wannier functions is 
the scaling of the original DFT calculations, which has 
been described elsewhere^i and we do not repeat now. 
However, this scaling has also an effect on the actual 
computation of WFs and we analyze this aspect in de- 
tail. 

As we described above, almost the whole amount of 
changes induced by the USPP description (relative to 
NCPP) in the the calculation of MLWFs is related to the 
implementation of Eq. H14|l . Details on how to compute 



III. APPLICATIONS 

In this Section we apply the above described formal- 
ism to some test cases, ranging from periodic crystals 
to isolated molecules. Precisely, we study fee Copper 
bulk, wurtzite AIN, and Watson-Crick DNA base pairs. 
We address a number of physical properties connected to 
Wannier functions: interpolation of the electronic struc- 
ture, calculation of dipole moments and spontaneous po- 
larization, analysis of the chemical bonding. All the cal- 
culations are performed with both norm-conserving and 
ultra-soft pseudopotentials. The numerical implications 
in using the USPP-TLA approach for the augmentation 
of overlaps are also discussed. 



A. Copper bulk 

We compute MLWFs for /cc-Copper, which has already 
been used as a test case in the literature^iSSiS of WFs. 
We adopt a 6 x 6 x 6 mesh of k-points to sample the Bril- 
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FIG. 1: (Color online) fee- Copper band structure: solid lines 
represent DFT Kohn-Sham bands (USPP) while dotted lines 
are results calculated using WFs. 



louin zone and compute six WFs corresponding to the 
lowest s ~ d manyfold. In the disentanghng procedure^ 
(used to get the optimal subspace for WF localization) 
we freeze the Bloch eigentates below the Fermi energy: 
this means that the subspace is constructed by these se- 
lected states plus a mixture of the states above the Fermi 
level. We adopt a kinetic energy cutoff for wavefuntions 
of 120 Ry (25 Ry) when using NCPP (USPP) and of 480 
Ry (200 Ry) for the density. In Fig. ^ we superimpose 
the band structure directly computed from a USPP-DFT 
calculation and that obtained from Wannier function in- 
terpolation^i^iii on the adopted 8x8x8 uniform k-point 
grid. These two sets of bands are almost superimposed 
below the Fermi energy and some slight differences arise 
only at higher energies. This is expected due to the choice 
of the energy window for the frozen states: while the 
eigenvalues for the k-points in the adopted regular mesh 
are the same as those from the DFT calculation by con- 
struction, it is not trivial that the band structure along 
a generic Brillouin zone line is well reproduced. This is 
indeed the case here, being a signature of proper localiza- 
tion of the computed WFs. The band structure interpo- 
lation obtained from NCPP is essentialy the same as the 
one in Fig. ^ and it is not reported. Some small descrep- 
ancies between NCPP and USPP interpolated bands are 
also present in the starting DFT calculations and are not 
of our interest in this context. 

In Tab.|l]we report a more detailed description of quan- 
tities related to WFs (spreads, real-space decay of the 
Hamiltonian matrix elements) in order to compare the 
NCPP and USPP approaches. A measure of the Hamil- 
tonian decay is defined as: 



d(R) 



l^^mn(R) 



1/2 



(15) 



where we defined iJ,„„(R) = (w;o,m I^^I^^r,™)- USPP re- 
sults appear to be in very good agreement with those 
related to NCPP, most of the differences being reason- 



TABLE I: fcc-Copper. WF spreads (Bohr^) and real -space 
decay of the Hamiltonian matrix elements (eV) for NCPP, 
USPP and USPP in the thermodynamic limit approximation 
(USPP- TLA). Q is the total spread, Qj, fin and Qqd are 
the invariant, diagonal and off-diagonal terms, according to 
Eq. iA3l . T is the (0.5 0.5 0.0) direct lattice vector. 





NCPP 


USPP 


USPP-TLA 


n 


22.434 


23.010 


18.886 




14.767 


15.629 


11.303 


^D+OD 


7.667 


7.381 


7.583 


d{ r) 


0.4876 


0.4795 


0.4779 


d(2T) 


0.0493 


0.0473 


0.0476 


d(3T) 


0.0203 


0.0207 


0.0205 



ably due to the pseudopotential generation and not to the 
WF computation. The average number of iterations to 
converge the disentanglement and the localization proce- 
dures are almost the same, as well as the singular spread 
values for each of the WFs. Figure |21 reports the spatial 
distribution of WFs from USPP calculations. As in pre- 
vious works)^*^ we find one interstitial s-like WF with 
the largest spread [Fig.|2a)], and five more localized d- 
like WFs [Fig.[2Ib-d)], correctly reproducing the physical 
s — d picture of Copper. 

As a last remark for the case of Copper, we analyze 
the effect of neglecting the e"''^ '" term in Eq. (|14|l . i.e. 
the USPP thermodynamic limit approximation (USPP- 
TLA). The theoretical background has already been dis- 
cussed in Sec. Ill CI here we focus on the numerical as- 
pects. In Tab. m the third column reports the resukts of 




(a) 



(b) 





(c) 



(d) 



FIG. 2: fee-Copper charge distribution for WFs computed 
using USPP. (a) interstitial WF, = 11.28 Bohr^; (b-d) d- 
charaeter WFs, = 1.86 Bohr^ Qc = 2.70 Bohr^, Qd = 1.54 
Bohr^. The two d-like missing WFs are strictly similar to (b) 
and (c) and are not reported. 
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the calculation performed within this approximation: it 
is evident that the numerical values of the USPP-TLA 
spreads deviate from the NCPP ones much more than 
the USPP do. On the contrary, the interpolated band 
structure and the real-space decay Hamiltonian matrix 
elements are definitely well-suited and comparable with 
those obtained in the full USPP treatment. Since the 
TLA is known to become exact in the thermodynamic 
limit, we expect it to work better when incresing the di- 
mension of the k-point mesh. We checked the behavior of 
the approximation with respect to different meshes but 
no convergence could be reached for grids ranging from 
4 X 4 X 4 to 10 X 10 X 10 k-points. 



B. Isolated molecules: DNA bases and base pairs 

Wannicr functions have been widely used to char- 
acterize the electrostatic properties of several molec- 
ular systems ranging from e.g. wates^SiSliiSi and 
small moleculesSiSi^iiS to large biomolecules, such as 
proteinSf^ nucleic acidsj^Si^ enzymes^ and ionic chan- 
nels.^*' In fact, the Wannier transformation allows one to 
partition the charge density into localized distributions of 
charges sitting on the so-called Wannier centers (f)„i^ 
In the case of isolated molecular systems, the dipole mo- 
ment is a well defined quantity given by p = pion + Pel 
with 

Pion ^ ^I'^i^ (16) 

I 

occ 

Pel = -2e^(r)„, 

n 

where pion and Pei are the ionic and electronic com- 
ponent respectively; e is the electron charge; the / sum- 
mation is over the ionic sites R/ and Zi is the valence 
charge of the /'^ atom, as defined by the correspond- 
ing psudopotential; the n summation is over the doubly 
occupied valence states. 

As a key test, we calculated the dipole moments of the 
four isolated DNA bases: Guanine (G), Cytosine (C), 
Adenine (A) and Thymine (T), and of the two Watson- 
Crick base pairs G-C and A-T, whose structures are re- 
ported in Fig. |3Ja). 

We simulated each system in a large (22 x 22 x 22) 
A'^ supercell, which allows us to avoid spurious interac- 
tions among neighbor replicas. For isolated systems, the 
uniform k-point grid reduces to the case of F-point only, 
and the connecting vectors b [Eq. |0J] correspond to the 
generators of the reciprocal lattice vectors. We expanded 
the electronic wavefuncions with the kinetic energy cutoff 
of 25 and 80 Ry, using USPP and NCPP respectively. 

We first optimized the atomic structure of the two base 
pairs G-C and A-T until forces on all atoms were lower 
than 0.03 eV/A, using ultra-soft pseudopotentials. Then, 



TABLE II: DNA bases: Dipole moments |p| (Debye) for 
isolated nucleobases (G, C, A, T) and for Watson-Crick 
base pairs (G-C, A-T). Present work results for both USPP 
and NCPP approaches are compared with previous quantum 
chemistry HF/6-31G** calculations-'^ and experimental data. 





USPP 


NCPP 


HF/6-31G" 


Exp. 


G 


7.1 


7.2 


7.1 


7.1" 


C 


6.7 


6.9 


7.1 


7.0'' 


A 


2.3 


2.3 


2.5 


2.5" 


T 


4.2 


4.4 


4.6 




G-C 


4.9 


4.9 


6.5 




A-T 


1.7 


1.8 


2.0 





"DeVoe and Tinoco (Ref.El) 
''Weber and Craven (Re f. 163) 
'^Kulakowski et al. (Ref. |6^ 



keeping atoms fixed in the relaxed geometry, we calcu- 
lated the electronic structure and the corresponding ML- 
WFs for both the single bases and the base pairs. We 
maintained the same geometries also for the correspond- 
ing NCPP calculations. In this case we have checked 
that forces on single atoms never exceeded the value of 
0.05eV/A. 

Our results for both sets of calculations are reported 
in Table m We note a very good internal agreement 
between the USPP and NCPP cases, as well as in com- 
parison with previous quantum chemistry HF/6-31G** 
calculationsSi and experimental results. The total spread 
O and its components f2/, Qqd (not shown) are also very 
similar in the two set of calculations, while the diagonal 
term is, by definition^ identically zero for isolated 
systems. 

Finally, we can take advantage from the calculated 
Wannier functions to further investigate the electronic 
distribution and the bonding pattern in the moleculeSi— 
For example, as shown in Fig. Ofb) for the case of the 
G-C base pair, we are able to characterize different kinds 
of bonds. The SB Wannier function represents a single a 
bond: It is centered in the middle point of the C-C bond. 
The two DB WFs describe instead a double N=C bond: 
This bond is partially polarized with a slight charge accu- 
mulation near the nitrogen atom. Finally, the LP WFs 
represents two electron lone pairs localized around the 
oxygen atom of the guanine molecule. Note also that 
the distribution of single/double bonds correctly reflects 
the theoretical one reported in Fig.|3fa). Let us remark 
here that our benchmark to assess the success of the 
newly developed USPP-WF methodology is its relative 
performance with respect to the NCPP-WF (already es- 
tablished) framework, namely the comparison between 
columns 2 and 3 in Table 11. 
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FIG. 3: (Color online) DNA bases: (a) Chemical scheme 
of the four DNA bases (G, C, A, T) assembled in the two 
Watson-Crick base pairs (GC and AT), (b) Isosurface plots of 
selected MLWFs for the GC base pair, representing a single 
C-C bond (SB), a double N=C bond (DB) and two oxygen 
lone pairs (LP). Dashed lines represent hydrogen bonds that 
bind the base pairs. 



C. AIN Wurtzite 

Here we move to the calculation of polarization in peri- 
odic systems. We focus on the case of aluminum nitride 
and compute the spontaneous polarization Pgp of the 
wurtzite phase. This is a particularly appealing test-case 
since a large debate exists in the literature and many 
results are present i2iSSiSLS2iS2ii2iIi Moreover, Nitrogen 
may be more easily described (in terms of PW kinetic- 
energy cutoff) with USPPs than with NCPPs, allowing 
for an advantageous application of the current formal- 
ism. Following Refs. ;9' 72] we evaluate the polarization 
(electronic and ionic contributions) for the wurtzite (WZ) 
phase taking the zinc-blend (ZB) structure as a reference. 
Our calculations adopt the GGA-PBE parametrization 
for the exchange-correlation functional, and uses cutoff 
energies of 60 Ry (25 Ry) for NCPP (USPP) for wave- 
functions and 240 Ry (200 Ry) for the density. We relax 
both the cell dimensions and the atomic positions of the 
WZ phase. The ZB reference is assumed to have ideal 
atomic positions and the cell is taken equal to the one 
computed for the WZ polytype (six bilayers are included 



TABLE IIL ^4^; electronic spontaneous polarization [C/m^] 
of wurtzite structure. The presented results (NCPP, USPP 
and USPP-TLA) are compared to the literature. As a refer- 
ence we also report the invariant (fli) and the total spread (Q) 
[Bohr^] given by the WF calculations for wurtzite. Results 
from Ref. are obtained using GGA and LDA (in parenthe- 
ses) respectively. 



Psp 


NCPP 


USPP 


USPP-TLA 


this work 


-0.095 


-0.094 


-0.094 


Ref. [9] 






-0.090 (-0.099) 


Ref. 




-0.120 




n 


80.895 


80.756 


80.271 




67.683 


67.857 


67.282 



in the cell) . These structural calculations have been per- 
formed using USPP and the obtained (lattice and ionic) 
parameters have been used also in the NCPP simulations. 
Using a 16 X 16 X 4 Monkhorst-Pack^^ mesh of k-points 
we obtained the a and c/a parameters of the exagonal 
lattice as a=3.1144 A and c/a=1.6109 (for the standard 
WZ cell including two bilayers, c/a=4.8326 for the actual 
cell we adopted). 

In Tab. IIIII we report the comparison of USPP and 
NCPP results in our calculations as well as other results 
from the literature. The values of spontaneous polar- 
ization computed using NCPPs are almost identical to 
those obtained with USPPs (Psp =-0.094 C/m^). The 
USPP-TLA behaves very accurately in this case and no 
difference can be found with respect to the full USPP 
calculation. Our results are also in very nice agreement 
with the USPP-TLA calculation by Bernardini et al.^ 
who used a Berry-phase formalism^ and found a value of 
Psp =-0.090 C/m^. We thus conclude that the USPP- 
TLA approximation performs well for the computation 
of the spontaneous polarization in nitrides, relative to a 
pure USPP treatment without the thermodynamic limit 
approximation. 

These results for the Psp in AIN wurtzite are also in 
agreement with __indirect experimental evidencies as re- 
ported in Refs. 68 69 7l]. While some earlier experimen- 
tal fits^'*''^''^ claim for much lower values of Psp (ranging 
from -0.040 to -0.060 C/m-2), later^LIL works explain 
this discrepancy as due the neglecting of bowing effects 
(non-linearity) of the Psp with respect to the composi- 
tion of the alloy Ala;Gai_2;N employed in the measure- 
ments. 



IV. CONCLUSIONS 

In this paper we presented an approach to calculate 
maximally localized Wannier functions in the ab ini- 
tio plane-wave ultrasoft pseudopotential scheme. Our 
methodology is formulated in the general framework of 
the PAW theory and recovers the USPP framework as 
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a special case. The main advantage using USPP is the 
weh known reduction of the computational effort in the 
evaluation of the electronic wavefunctions at the DFT 
level. This leads to a consequent reduction of computa- 
tional load also in the calculation of MLWFs. We demon- 
strated that the extension to the USPP case does not 
introduce further approximations in the computation of 
the MLWFs with respect to the NCPP case. Further- 
more, the reformulation within the PAW scheme allows 
us to interface the computation of MLWFs to other pop- 
ular approaches for the electronic structure calculation. 
Finally, our method is formulated in the case of a uni- 
form k-point mesh for the sampling of the Brillouin Zone, 
generalizing previous attempts based on P-only calcula- 
tions. This allows us to treat periodic solid-state systems 
(such as crystals, surfaces and interfaces), which require 
a full description of the BZ, as well as molecular, finite 
or amorphous systems which are well described with the 
F-only representation. 

As a first illustration of the capability of this method- 
ology, we presented the calculation of the MLWFs for a 
few selected test cases, easily referable to well established 
theoretical and experimental results. For each selected 
system we also compared the results for both USPP and 
NCPP calculations, underlying a very good agreement 
between the two cases. 

The reduction of the computational cost resulting from 
USPP calculations opens the way to the exploitation of 
the MLWFS as a powerful tool to analyze the electronic 
structure of larger and more realistic nanoscale systems, 
in particular for transport in nano-j unctions . ^ ^ i ^ ^ 
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APPENDIX A: MAIN FORMULAS FOR 
MAXIMALLY LOCALIZED WANNIER 
FUNCTIONS 

For sake of completeness we report here the main 
relations^ entering the expression (|3Jl for the spread func- 
tional Q [U] and its first derivative wrt the unitary trans- 
formation U [Eq. These are all the quantities in- 
volved in the minimization procedure for the calculation 
of maximally localized Wannier functions. The expecta- 



tion values of the position operator are: 

= -^Y.'"bhlm\nM^;^, (Al) 

k,b 

k,b 

where b vectors connect nearest-neighbor k-points and 
Wb are their weights according to Appendix B of Ref. ■ 
The spread functional can be divided into three terms, 
the invariant (I), the diagonal (D) and the off-diagonal 
(OD) components:'^ 

n[u] ^ ni + noiu] + noD[u], (A3) 

Their definitions are, respectively: 




= ]^ E E (lmlnM„^^ + b . r„)' , (A5) 

k,b n 

APPENDIX B: DETAILED EXPRESSIONS FOR 
THE CALCULATION OF AND Q»j(b) 

We report the explicit expression for the calculation 
of the reciprocal space representation of the PAW/USPP 
projectors and the Fourier transform of the augmentation 
densities: 

Qi.ih)^ J dve-^'^-^Ql^iv-n). (Bl) 

These tasks are also required in the evaluation 6.17. of the 
density in reciprocal space and are therefore performed 
in standard plane-waves DFT codes. The index j,j in 
I3i(v) and Qiji^) stand for radial and angular numbers, 
6.17. niliUii. Projectors and augmentation densities are 
explicitly written as a product of a radial part times (real) 
spherical harmonics: 

A(r) = i?„.(r)l^™'(f) (B2) 
Q,,{r) = ff„.„,(r)l^™'(r)1^7(f). (B3) 

First we focus on the expression for /3 projectors. Func- 
tions of the form /(r) = R{r) Y^^f) have a known semi- 
analytical Fourier transform which is given by: 

/>oo 

/(k) =47r(-i)'y™(fc) / R{r) Ji{kr)dr (B4) 
Jo 

where Ji{x) is the spherical Bessel function of order 
The problem for (3 projectors is therefore directly solved 
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once we add the structure factors accounting for the acu- 
tal positfons of the fon: 



X 



X / R,,^{r)Ji^{Gr)dr. (B5) 
/o 



Moving to the case of the augmentation densities, we 
note that the product of two spherical harmonics can 
be expressed as a hnear combination of single spherical 
harmonics using Clebsch-Gordan coefficients: 



L=\li-lj\ M=-L 



(B6) 



This allows to follow the same strategy as before also 
for Eq. lUll). Putting Eqs. (|B6IB4(I together, the final 
expression for Qjj{h) reads: 



L M 

xYi'ib) / Qn^nM J L{hr)dr. (B7) 
Jo 



where M, L indexes run as in Eq. ((B 
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